The double-funnel energy landscape of the 38-atom Lennard-Jones cluster 
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The 38-atom Lennard-Jones cluster has a paradigmatic double-funnel energy landscape. One 
funnel ends in the global minimum, a face-centred-cubic (fee) truncated octahedron. At the bottom 
of the other funnel is the second lowest energy minimum which is an incomplete Mackay icosahedron. 
We characterize the energy landscape in two ways. Firstly, from a large sample of minima and 
transition states we construct a disconnectivity tree showing which minima are connected below 
certain energy thresholds. Secondly we compute the free energy as a function of a bond-order 
parameter. The free energy profile has two minima, one which corresponds to the fee funnel and the 
other which at low temperature corresponds to the icosahedral funnel and at higher temperatures 
to the liquid-like state. These two approaches show that the greater width of the icosahedral funnel, 
and the greater structural similarity between the icosahedral structures and those associated with 
the liquid-like state, are the cause of the smaller free energy barrier for entering the icosahedral 
funnel from the liquid-like state and therefore of the cluster's preferential entry into this funnel 
on relaxation down the energy landscape. Furthermore, the large free energy barrier between the 
fee and icosahedral funnels, which is energetic in origin, causes the cluster to be trapped in one of 
the funnels at low temperature. These results explain in detail the link between the double-funnel 
energy landscape and the difficulty of global optimization for this cluster. 



I. INTRODUCTION 

Understanding the relationship between the potential 
energy surface (PES), or energy landscape, and the dy- 
namics of a complex system is a major research effort in 
the chemical physics community. For example, much the- 
oretical work has attempted to find the features of the en- 
ergy landscape which differentiate those model polypep- 
tides that are able to fold rapidly to a unique native striic=. 
ture from those that get stuck in misfolded states.ETEl 
Similarly, the answers to a whole host of questions about 
glasses lie in the energy landscape. Why are glasses un- 
able to reach the crystalline state? What is the cause 
of the differences between 'strong' and 'fragile' liquids? 
What processes, at a microscopic level, are responsible 
for a and [3 relaxation? 

A key concept that has arisen within the protein fold- 
ing community is that of a funnel consisting of a set of 
downhill pathways that converge on a single low-energy 
minimum.Qu It has been suggested that the PES's of pro- 
teins are characterized by a single deep funnel and that 
this feature underlies their ability to fold to their native 
state. Indeed it is easy to design model single-funnel 
PES's that result in efficient relaxation to the glp.Hal 
minimum, despite very large configurational spaces.ErtL£l 
By contrast, polypeptides that misfold are expected to 



have other funnels that can act as traps. Attempts 
have been made to characterize the energy landscape of 
proteins in these terms through, for example, mapping 
the connections -between compact states, nn 1 j disconnec- 
tivity graphsjjjjwp monotonic sequenceaij and free en- 
ergy profiles .Eillia However, for the most part these stud- 
ies have been limited cither to simplified lattice models, 
where the most natural elementary division of the .EES 
into basins of attraction surrounding local minimatZI is 
problematic, or to short polypeptides if more realistic 
models are used. 

An intuitive picture of the energy landscape of glasses 
has also been proposed in which the crystal corresponds 
to a vesy-jnarrow funnel which is inaccessible from the 
liquid Jl§ll2l Most of the PES is dominated by rugged re- 
gions where there are many funnels leading to different 
amorphous structures. Therefore, the structure that the 
system relaxes to depends upon its thermal history. In 
this picture the different relaxation processes result from 
the hierarchy of barriers in the amorphous regions of the 
PESO Although this picture is appealing, only recently 
has progress been made in relating the-details of glassy 
behaviour to the features of the PES .£21123 This task is 
hampered by the complexity of the PES's for these sys- 
tems; the number of minima is huge and characterization 
of the PES is hampered by the slow relaxation times. 
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In view of the complex energy landscapes of model 
proteins and glasses, it is also useful to have systems 
for which it is easier to understand the relationship be- 
tween the PES and the dynamics. Clusters can provide 
just such an alternative perspective. For example, the 
complexity of the PES (in terms of the numbe r o f .min- 
ima) can be controlled through the cluster sizecaEJ and 
the choice .of. potential parameters, such as the range of 
attraction.EjuJ 

Clusters where the atoms interact through the 
Lennard-Jones (LJ) potential — 



E 



i<j 



(1) 



where e is the pair well depth and 2 1 / 6 <r is the equilibrium 
pair separation — provide a particularly useful model sys- 
tem, because their structure, thermodynamics and dy- 
namics have been much studied. For small LJ clus- 
ters a complete enumeration of the minima and tran- 
sition states-allows a detailed view of the dynamics to be 
obtamedBMI 

At larger sizes, well-chosen examples allow one to con- 
sider particular paradigmatic types of energy landscape. 
For example, the PES of LJ55 has a single deea funnel 
which leads down to the Mackay icosahedror£3 global 
minimum, and LJ38, the cluster which we study here, has 
a double-funnel landscape. The LJ38 global miBJ«ium is 
a face-centred-cubic (fee) truncated octahedronEffill (Fig- 
ure |l|a) and the second lowest-energy minimum is an in- 
complete Mackay icosahedronEa (Figure [l]b). These two 
minima lie at the bottom of separate funnels on the PES. 
The thermodynamics of this cluster have recently been 
characterizedHo At low temperature the truncated oc- 
tahedron has the lowest free energy but, because the en- 
tropy associated with the icosahedral funnel is larger, its 
free energy becomes lower than that for the fee funnel at 
about two thirds of the melting temperature. Therefore, 
a transition takes place between the two states which is 
the finite-size equivalent of a solid-solid phase transition. 




FIG. 1. (a) The LJ38 global minimum, an fee truncated 
octahedron (E = -173.928427e; point group O h ). (b) Sec- 
ond lowest energy minimum of LJ38 (E — — 173.252378e; 
point group Csv)- (c) Third lowest energy minimum 
(E = — 173.134317e; point group C s ). The structures in 
(b) and (c) are both incomplete Mackay icosahedra. 



This thermodynamic transition affects the dynam- 
ics. Relaxation down the PES from the liquid-like 
state almost invariably leads into the icosahedral fun- 
nel. This is partly because, near to melting, icosahe- 
dral structures have a lower free energy than fee struc- 
tures. However, entry into the icosahedral funnel also 
seems to be dynamically favoured, perhaps because of 
the greater structural similarity between the icosahedral 
and liquid-likejKtniatures — both have some polytetrahe- 
dral character .QOO One aim of this paper is to charac- 
terize in detail the reasons for the greater accessibility of 
the icosahedral funnel. 

Furthermore, once the cluster enters the icosahedral 
funnel it becomes trapped there, even when the free en- 
ergy of the fee funnel is lower. There is a large free energy 
barrier between the two funnels which prevents the clus- 
ter passing between them, but the nature and size of this 
barrier have not yet been probed. 

These features make global optimization of this sys- 
tem much more difficult than for most small LJ clusters; 
most have global minima based on the Mackay icosahe- 
dra with no other competitive morphologies x3 Indeed, 
the LJ38 global minimuift,was initially discovered on the 
basis of physical insight .E3 Although the truncated octa- 
hedron has since beaa . fauj id-by a number of global op- 
timization methods J^'EBEJe3 most of these techniques 
examine not the usual LJ energy landscape, but a land- 
scape that has been transformed with the aim of making 
global optimization easier. The basis for the success of 
one of these methods lies in the significant changes to the 
thermodynamics and dynamics that the transformation 
causes Ma 

We use two tools to characterize the double-funnel to- 
pography of the LJ38 PES: in section [n| we use a discon- 
nectivity graph, and in section |l| we examine the free 
energy landscape. Some of the results presented here 
have previously appeared in a short communication.EU 



II. DISCONNECTIVITY GRAPH 

To examine the topography of the PES, we need to 
locate its minima and the network of transition states 
and pathways that connect them. A transition state is 
a stationary point on the PES where the Hessian matrix 
has exactly one negative eigenvalue. Transition states 
can be found efficiently using eigenvector-followingJUo 
in which the energy is maximized along one direction 
and simultaneously minimized in all others. The minima 
connected to a transition state are defined by the end- 
points of the steepest-descent paths commencing parallel 
and antiparallel to the transition vector (the eigenvec- 
tor whose eigenvalue is negative) . For this calculation we 
employ a method that uses analytic second derivatives .e3 

The number of locally stable structures that LJ38 can 
adopt is too large for it to be desirable or even possible to 
catalogue them all. However, here we are primarily inter- 
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ested in the energetically low-lying regions of the PES as- 
sociated with the two funnels. In recent years, a number 
of similar approaches for systematically exploring-a-EES. 
by hopping between wells have been developedpj'EJCj 
and these are easily adapted to produce an algorithm that 
explores low-energy regions of the PES preferentially. In 
our scheme, we commence at a known low-lying minimum 
and proceed as follows. 

1. Search for a transition state along the Hessian 
eigenvector with the smallest non-zero eigenvalue. 

2. Find the steepest-descent pathway through the 
transition state and the two minima it connects, 
as described above. 

3. There are various possible outcomes from step 

(a) In most cases, one of the connected minima is 
the minimum from which the transition state 
search was initiated. If this is the case, and the 
other minimum is lower in energy, we move to 
it. 

(b) If the original minimum is one of the con- 
nected ones, but the other is higher in energy, 
the move is rejected. 

(c) Sometimes the transition state is not con- 
nected to the minimum from which the search 
started. If neither minima has been found pre- 
viously, the pathway is then isolated from the 
rest of the database. Since we want to explore 
patterns of connectivity in the low-energy re- 
gions of the PES, we discard the transition 
state and both minima under these circum- 
stances. Such searches can be repeated later, 
when the database has grown and a connec- 
tion may be found. 

(d) If the original minimum is not connected, but 
one or both minima have already been visited, 
the pathway is recorded, but we remain at the 
original minimum. 

4. The procedure continues from step |], searching in 
both directions along eigenvectors with successively 
higher eigenvalues from the modified position. 

5. When a specified number, n CVl of eigenvectors of a 
minimum have been searched for transition states, 
the position jumps to the lowest-energy known min- 
imum for which fewer than n cv eigenvectors have 
been searched. 

By only accepting downhill moves, this algorithm pre- 
vents the search becoming lost in the manifold of liquid- 
like minima. In the present work, we chose n cv = 10, 
allowing up to 20 transition state searches from each 
minimum. In general, searches along eigenvectors with 
low eigenvalues are more likely to converge to transition 



states in a reasonable number of iterations. One can ob- 
tain an impression of how thoroughly the low-energy re- 
gions of the PES have been explored by monitoring how 
many of the lowest-energy minima are displaced as the 
search proceeds. Having collected 3000 minima, the next 
1000 displaced only 6 of the previous lowest 200, and the 
next 2000 displaced only a further 7. 

In the course of the search, several multiple-step paths 
between the low-lying Oh and C$ v minima emerged. The 
highest point on the lowest-energy path that we found 
was a transition state with energy — 169.709e (4.219e 
above the global minimum); the path has an integrated 
length of 18.517cr. This demonstrates the efficiency of 
well-hopping techniques over more conventional methods 
for exploring a PES, such as molecular dynamics (MD): 
to restrict MD to regions of low potential energy, the to- 
tal energy of the simulation would have to be low, and the 
trajectory would waste time undergoing intrawell oscil- 
lations, but at energies high enough to allow interfunnel 
passage at a reasonable rate, the trajectory could escape 
into the numerous liquid-like structures. 

One way to analyse a database of minima an d t ran -. 
sition states is in terms of monotonic sequences ,EJc3'E2l 
i.e. connected sequences of minima with monotonically 
decreasing energy, which terminate at a particular min- 
imum. A funnel can then be defined as containing all 
minima which lie on monotonic sequences to the lowest 
common minimum. The primary funnel terminates at 
the global minimum, and is separated from secondary 
funnels by primary divides — higher-lying minima lying 
at the boundary between funnels. Classifying regions of 
the PES in this way is useful because motion between 
minima in a given funnel is likely to rQCCiir on a shorter 
time scale than flow between funnels. Cj'Eil 

One can also gain insight into dynamical features of a 
system by examining the energy profile of the monotonic 
sequences. If the decrease in energy between minima in 
the sequence (i.e. the energy gradient towards the funnel 
bottom) is large compared with the intervening barri- 
ers, so that the profile looks like a staircase, relaxation 
towards the funnel bottom is likely to be easy. This ap- 
proach has been used to explain the structure-seeking 
properties of a KC1 cluster .l3 In contrast, energy profiles 
that have a shallow gradient and high barriers are more 
likely to produce glass- formers. 

An informative way to visualize the PES and-^feveal 
funnel structure is to plot a disconnectivity graph.tHl Such 
graphs have been used to obtain—insight into the en-, 
ergy landscape of polypeptidesJi3t3 Cgo^l and water 
sodium chlorideEj and MorseEj clusters. At a given to- 
tal energy, E, minima can be grouped into disjoint sets, 
called basins, whose members are mutually accessible at 
that energy. In other words, each pair of minima in a 
basin are connected directly or through other minima by 
a path whose energy never exceeds E, but would require 
more energy to reach a minimum in another basin. At 
low energy there is just one basin — that containing the 
global minimum. At successively higher energies, more 
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basins come into play as new minima are reached. At still 
higher energies, the basins coalesce as higher barriers are 
overcome, until finally there is just one basin containing 
all the minima (provided there are no infinite barriers). 

In a disconnectivity graph, the basin analysis is per- 
formed at a series of energies, which are plotted on a 
vertical axis. At each energy, a basin is represented by 
a node, and lines join nodes in one level to their daugh- 
ter nodes in the level below. The horizontal position of 
a node has no significance, and is chosen for clarity. A 
funnel appears as a tall stem with branches sprouting 
directly from it at a series of levels, indicating the pro- 
gressive exclusion of minima as the energy is decreased. 
In the language of Ref. |Il], this pattern resembles a palm 
tree. Significant side-branching would indicate additional 
features of the PES and leads to different patterns .ell 
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FIG. 2. Disconnectivity graph for LJ38 using a sample of 
6000 minima and 8633 transition states. Only branches lead- 
ing to the 150 lowest-energy minima are shown, but numbers 
attached to nodes indicate the number of minima they rep- 
resent. The branches terminating at the three lowest-lying 
minima (see Figure ^) are labelled by their point groups. 
The energy scale is in units of e. 

Figure |^ shows the disconnectivity graph for LJ38 using 
a sample of 6000 minima and 8633 transition states, and 
a level spacing of 0.5e. For clarity, only branches leading 
to the 150 lowest-energy minima are shown, but the num- 
ber of minima that would be represented by some of the 
larger nodes are indicated. The large funnel associated 



with the C§ v minimum (placed centrally) is immediately 
visible. The lowest node connecting it to the funnel of the 
fee global minimum lies at — 169. 5e, so the 446 minima 
in the node below can be considered as belonging to the 
icosahedral funnel. In contrast, the funnel of the global 
minimum contains only 47 of the minima in our sample — 
nearly an order of magnitude fewer. Although only 150 
branches are shown in the figure, one can see from the 
rapidly increasing density of branch-ends as the energy 
rises past — 171. 5e, that the number of states available 
to the cluster increases dramatically with energy. This 
increase signifies the onset of the liquid-like part of the 
PES, and the disconnectivity graph shows that the sys- 
tem must enter this region of configuration space in order 
to pass between the icosahedral and fee funnels. 

The graph also gives an impression of the "shape" of 
the two funnels. The branches in the fee funnel are gen- 
erally longer than those in the icosahedral funnel, in- 
dicating higher barriers. Also, the global minimum is 
considerably lower in energy than the rest of the minima 
in the fee funnel because it has a complete outer shell. 
In contrast to the fee funnel, and to sizes at, which a 
complete Mackay icosahedron can be formed,c3 the bot- 
tom of the icosahedral minimum is not dominated by 
a single minimum. Rather there are many low-energy 
icosahedral minima associated with different ways of ar- 
ranging the atoms in the incomplete surface layer. It is 
also noticeable that the barriers about the low-lying C s 
minimum (Figure Qc) are lower than that for the C§ v 
minimum indicating that rearrangement of the vacancy 
associated with the missing vertex atom is much easier in 
the C s minimum. These features probably explain—why 
optimization methods found the C a minimum first ;E3 the 
C$ v minimum was only discovered relatively recentlyEj 

The overall picture, therefore, is of a narrow, deep 
and somewhat rougher funnel containing the global min- 
imum, and a broader, much more voluminous funnel as- 
sociated with the low-lying icosahedral minimum. The 
"rims" of both funnels lie in the liquid-like regions of the 
PES. The greater width of the icosahedral funnel helps to 
explain why the cluster enters this region of configuration 
space in the vast majority of annealing simulations. This 
effect of the funnel width has bjsen previously observed 
for a model double- funnel PESJa 

Our sample of minima is clearly only a tiny fraction 
of the astronomical number available to the system, so 
we need to consider the possible effects of incomplete- 
ness for the disconnectivity graph. As discussed above, 
from about half way through the search, very few new 
structures had lower energy than any of the existing low- 
est 200. This provides good evidence that the 150 min- 
ima actually represented in Figure || really are the low- 
est. The number of minima represented by higher-energy 
nodes (for which branches have not been shown) would 
certainly increase if the search were allowed to proceed 
for longer. 

The incompleteness of the transition state sample is 
harder to gauge, but has important consequences for the 
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disconnectivity graph. Given an incomplete sample of 
minima, the graph depends only on transition states that 
interconnect minima within the sample. Furthermore, 
two minima may be connected by more than one tran- 
sition state, but only the lowest matters for the graph 
because it determines the energy at which the minima 
become mutually accessible. When a new connectivity 
is discovered, the pattern of nodes and lines may change 
significantly. When a lower transition state between two 
previously connected minima is discovered, branching 
moves down the graph to lower nodes. In the present 
work, up to 20 transition state searches were allowed 
from each minimum, with many of the low-energy min- 
ima reaching this limit. A total of 25 403 transition state 
searches were performed, 21 515 of which converged in 
a reasonable number of optimization steps. Since only 
8 633 transition states were found, most of them must 
have occurred several times, suggesting that we have an 
adequate sample, especially as far as the low-energy min- 
ima are concerned. 

In summary, therefore, we can be quite confident that 
the disconnectivity graph in Figure is an accurate rep- 
resentation of the low-energy regions of the PES. This 
is because it is based on a search that is strongly biased 
towards low-energy structures, and because the sample 
of minima and transition states is far larger than the 
number of branches actually included on the graph. 



III. FREE ENERGY LANDSCAPE 

An alternative way of characterizing the double-funnel 
topography of the LJ38 PES is to compute the free energy 
as a function of a suitable order parameter. The two fun- 
nels would be expected to give rise to minima in the free 
energy which are separated by a barrier. Bond-order pa- 
rameters, which were initially introduced by Steinhardt 
et al.,E3 have been used to characterize the frefi-.epergy 
barrier for the nucleation of a crystal from a meltEiTEJ be- 
cause they can differentiate between the fee crystal and 
the liquid. By calculating the bond-order parameters, 
Qi, Qg, W4 and Wq, for our sample of LJ38 minima we 
were able to assess whether they might also be used to 
differentiate the fee and icosahedral funnels. Both Q4 
and Qq appeared suitable and we chose to investigate Q4 
further. 

The definition of the order parameter, Qi, 



is 



\ m— — l / 



where 



(2) 



(3) 



nearest- neighbour criterion, ro) in the cluster, Yi m {Q,<j>) 
is a spherical harmonic and 0y and fyj are the polar and 
azimuthal angles of an interatomic vector with respect to 
an arbitrary coordinate frame (Qi is independent of the 
choice for this coordinate frame). We use r — 1.391cr. 
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where the sum is over all the Nb 'bonds' (pairs of atoms 
which have a pair separation, r^, which is less than the 
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FIG. 3. Qi and potential energy values for a set of 7000 
LJ38 minima. The points associated with the three low- 
est-energy minima are labelled by their point group. 



Figure || shows that the Q4 values of the two lowest 
energy minima are well-separated. However, the large 
number of minima associated with the liquid-like state 
have values only slightly greater than those for the icosa- 
hedral structures. Therefore, Q4 is a good order param- 
eter for distinguishing fee structures but not for differen- 
tiating the icosahedral and liquid- like structures .Ell The 
similar values of Q4 for the icosahedral and liquid-like 
minima reflects structural similarities — both have signif- 
icant polytetrahedral character .EjEj 

In Figure ^ we show the properties of the lowest-energy 
pathway between the two lowest-energy minima that we 
found in section ||. The value of Q4 rapidly decreases 
as the cluster leaves the fee funnel and enters the liquid- 
like region of configuration space. A number of rear- 
rangements take place between liquid-like minima before 
the cluster then enters the icosahedral funnel. From the 
pathway we can obtain upper bounds to the free energy 
barriers between the two funnels at zero temperature: the 
T = free energy barriers must be less than 4.22e and 
3.54e for fee to icosahedral and icosahedral to fee tran- 
sitions, respectively. These values are upper bounds be- 
cause there may be lower-energy pathways that we were 
unable to find. Also the free energy difference between 
the two minima at T — is simply the energy difference, 
0.68e. 

To confirm that Q4 is not only able to distinguish fee 
and icosahedral minima, but also more general configu- 
rations from within the two funnels, we performed two 
series of Monte Carlo simulations of increasing temper- 
ature that started from the two lowest energy minima. 
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The probability distributions of Q4 for the two runs are 
well separated and do not overlap until until the cluster 
starts to melt at T ~ 0.17efc _1 (Figure ||). Below the 
melting temperature the clusters remain in the funnel in 
which the simulations were started even when the other 
funnel has a lower free energy. Hence there is a free en- 
ergy barrier between the two funnels which is significantly 
larger than the thermal energy. 
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FIG. 4. (a) Q4 as a function of the distance along the 
reaction pathway (starting at the global minimum) and (b) 
Qa against potential energy for the lowest energy pathway 
between the two lowest energy LJ38 minima. The diamonds 
are for the stationary points (minima and transition states) 
on the pathway. 



In the canonical ensemble the Landau free energy is 
related to the equilibrium probability distribution of the 
order parameter by 



A L {Qi) = A-kTlogp^iQi), 



(4) 



where A is the Helmholtz free energy. However, con- 
ventional simulations are unable to provide equilibrium 
probability distributions for Q4 because, as Figure |^ il- 
lustrates, the cluster is unable to pass over the free en- 
ergy barrier between the two funnels. To overcome this 



difficulty we use umbrella samplingB In this method 
configurations are not sampled with a Boltzmann dis- 
tribution but with the distribution exp(— (3E + W(Qa)) 
where W(Qn) is a biasing potential, the aim of which is 
to make configurations near the top of the free energy 
barrier more likely to be sampled. The canonical proba- 
bility distribution is then obtained from the probability 
distribution for the biased run, p mu iti(Q4), by 



PcaniQi) =_Pmulti(<24)exp(-PF(<2 4 )). 



(•5) 



We wish to choose W such that Pmuiti(Q4) is approxi- 
mately constant over thepwiiole range of Q 4 (the so-called 
multicanonical approachP^I 60 !) . However, this only occurs 
when W{Qi) ~ A^Q^/kT and so we have to construct 
W iteratively from the. results of a number of short pre- 
liminary simulations .Ell 
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FIG. 5. Q4 for two series of canonical Monte Carlo runs 
of increasing temperature. The initial configurations were 
the fee global minimum (solid line) and the lowest-energy 
icosahedral minimum (dashed line). Each point is the aver- 
age value in a 2 x 10 6 cycle Monte Carlo run and each run 
was initiated from the final geometry in the previous lower 
temperature run. The error bars represent the standard de- 
viation of the Q4 probability distributions. 



We were able to compute Al{Qa) successfully for 
T > l.Sefc . However, at lower temperatures, even with 
a reasonable biasing distribution, the rate at which the 
system passed between the two funnels was too low for 
accurate free energies to be obtained. At the top of the 
free energy barrier in this temperature range the contri- 
bution of states which mediate transitions between the 
two funnels is small (instead the contribution of other 
lower-energy states dominates) , and so even at the top of 
the barrier a transition to the other funnel is an activated 
process. To overcome this difficulty would require either 
unfeasibly long simulations or a better order parameter 
for which the contribution of irrelevant states to the free 
energy barrier region is less. However, the results we ob- 
tain are sufficient to give us a good picture of the free 
energy landscape of L J38 . 
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In our simulations we collected the values of Q4 and 
E c (the configurational, or potential, energy) into a two- 
dimensional histogram. This approach allows us to ob- 
tain the two-dimensional free energy surface, Al(Q±, E c ) 
(Figure^), as well as the free energy profile, Al^Q^) (Fig- 
ure 0a). Furthermore, it allows us to decompose Al(Qh) 
into its energetic and entropic components, 

A L {Q i )=E c<L {Q i )-TS L {Q A ), (6) 

because we can obtain E c ,l(Q4) from our two- 
dimensional probability distribution using 

Ec,L(Qi)= I p can (Qi,E c )E c dE c . (7) 

Moreover, we can apply the histogram methodEl to cal- 
culate results for temperatures other than those at which 
we performed simulations. As 

p(Q 4 ,£ c ;/3) = tt c (Q 4 ,£ c )exp(-/3£ c )/Z(/3), (8) 

where n c (Qi, E c ) is related to the configurational den- 
sity of states by Sl c (E c ) = J fl c (Q4, E c )dQ± and Z{0) is 
the partition function, it follows that 

p(Qi, E c ; /3') oc P (Q 4 , E c ; /3) e X p(-E c (/3' - /?)), (9) 

where (3 is the reciprocal temperature of the original sim- 
ulation and [3 1 is the reciprocal temperature to which the 
results have been extrapolated. This method, though, 
has to be applied with a certain amount of caution since 
the extrapolation becomes increasingly sensitive to sta- 
tistical errors at the edge of the p{Qa, E c ; f3) distribution 
(Figure |^) as the difference in temperature increasesO 

At T = 0.15 e/c" 1 the free energy profile has two main 
minima (the one at Q4 — 0.015 corresponding to the 
icosahedral funnel and the one at Q4 = 0.186 correspond- 
ing to the fee funnel) which are separated by a barrier 
(Figure ^1). Around the fee free energy minimum there 
are a number of oscillations in Al^Q^). These are a result 
of the discontinuities that occur in the order parameter 
when the value of an interatomic distance passes through 
ro and do not indicate that there are small free energy 
barriers between different fee structures. These disconti- 
nuities can also been seen in the pathway in Figure || at 
large values of Q4. 

At T = 0.15 efc -1 the icosahedral funnel is lower in free 
energy because it has a larger entropy (Figure 0c and ||a) 
due to the larger number of icosahedral minima-CEjgure 
H) and their lower mean vibrational frequencyEjta The 
free energy barrier is large with respect to the thermal 
energy (11.5 and 9.69 kT with respect to the free energy 
minima) and relative to kT increases rapidly as the tem- 
perature decreases (Figure [)]b) . The size of the barrier ex- 
plains why, below the melting temperature, simulations 
are trapped in one well or the other. 
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FIG. 6. Contour plots of Al(Qa, E c ) at temperatures of 
(a) 0.15efc _1 , (b) 0.18efc _1 and (c) 0.21 efc- 1 . The contours 
are spaced 1 kT apart. The free energy zero is the bottom of 
the icosahedral/liquid-like free energy minimum. Near to the 
free energy transition state some of the contours have been 
labelled by the value of the free energy in kT. No points 
were sampled in the regions free of contours. The results are 
from MC runs of (a) 250 x 10 6 , (b) 20 x 10 6 and (c) 15 x 10 6 
cycles. 
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FIG. 7. (a) A L (Q 4 ) (b) E C , L (Q 4 ) and (c) S L {Qi) at three 
different temperatures, as labelled (in units of efc -1 ). In each 
case, we have chosen the low Q4 free energy minimum as the 
zero of Al and the lowest value of Sl in the fee funnel as 
the zero of Sl ■ 



FIG. 8. The decomposition of the free energy profiles 
Al(Q4) (solid line) into their energetic {E c ,l(Qa) [dashed 
line]) and entropic (— T ] Sl{Qi) [dotted line]) contributions 
at temperatures of (a) 0.15 efc -1 , (b) 0.18 efc -1 and (c) 
0.21 efc -1 . The zeros of all three quantities have been set 
to occur at the position of the low Q4 free energy minimum. 
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(a) 



0.1 0.12 0.14 0.16 0.18 0.2 0.22 
temperature /ek"' 




0.14 0.16 0.18 
temperature /£k"' 



(c) 




0.05 0.1 0.15 0.2 

FIG. 9. (a) and (b) The height of the free energy bar- 
rier relative to the fee (f) and low Q4 (i) free energy minima 
and the free energy difference (dotted line) between the two 
free energy minima as a function of temperature. The points 
are derived from the simulation profiles and the lines result 
from using a histogram method, (c) Two-dimensional con- 
tour plot showing dependence of the free energy landscape, 
Al(Qa, T) on temperature. At each temperature we have set 
the zero of the free energy to that of the lower free energy 
minimum. The units of free energy are e in (a) and (c), and 
kT in (b). The contours in (c) are at a spacing of 0.2 e 



Examining Figure 0a shows that at T = 0.15 ek^ 1 the 
free energy barrier is energetic in origin, and that the 
larger entropy of the intermediate states acts to reduce 
the magnitude of the barrier. As the temperature de- 
creases, this entropic contribution decreases and so the 
barrier (in absolute terms) increases until it reaches its 
purely energetic value at zero temperature. The his- 
togram approach allows an estimate of the position of 
the fee to icosahedral transition. It predicts that the 
free energy difference between the two funnels is zero at 
T = 0.118 ek- 1 (Figure |) which is in good agreement 
with the thermodynamic results obtained using the su- 
perposition method.Eil 

As the temperature increases from T = 0.15 e/c -1 the 
liquid-like state makes an increasing contribution to the 
low Q4 free energy minimum. For example, the value 
of Q4 at the minimum gradually increases (Figures 0a 
and 0c) , reflecting the slightly larger values of Q4 for the 
liquid-like state (Figure 0) and the minimum becomes 
broader. At T = 0.18 ek -1 in a canonical simulation 
the cluster passes back and forth between the liquid-like 
and icosahedral states in time. This dynamical coexis- 
tence is reflected in the low Q4 free energy minimum of 
Al(Q4, E c ) — comparison of Figure [|d to Figures 0a and 
0c shows that it is a superposition of two states. 

At T = 0.21 ek -1 the low Q4 f ree energy minimum is 
solely due to the liquid-like state and the free energy land- 
scape is dominated by the free energy difference between 
the fee and liquid-like structures which results from the 
much larger entropy of the liquid (Figure 0c). The fee free 
energy minimum is now shallow and the flatness of the 
free energy landscape for Q4 > 1.0 is a result of the com- 
pensation of the energy and entropic components (Figure 
0c). The histogram approach predicts that the fee free 
energy minimum finally disappears at T ~ Q.2ibek~ l 
(Figure|). 

The free energy barrier to pass from the low Q4 free 
energy minimum to the fee funnel has an interesting tem- 
perature dependence (Figure |). It has a minimum at a 
temperature close to the melting transition. Below this 
temperature the barrier increases because of the decreas- 
ing effect of the entropy of intermediate states in reducing 
the energetic barrier between the two free energy min- 
ima. Above this temperature the barrier increases be- 
cause of the increasing free energy difference between the 
two minima. As a canonical simulation is most likely to 
enter the fee funnel when the free energy barrier relative 
to the thermal energy is at its smallest, the optimum 
temperature for reaching the basin of attraction of the 
fee global minimum is T ~ 0.19efc _1 (Figure 0b). This 
result is somewhat counter-intuitive because one would 
usually expect the optimum temperature to be when the 
equilibrium probability of being in the basin of attrac- 
tion of the global minimum is higher: Pq — 0.004 at 
T = 0.19 ek^M Figure |l0| does indeed confirm that a 
simulation at this temperature can enter the fee funnel, 
albeit rarely. The cluster enters the fee funnel once in 
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1.5 x 10 7 MC cycles and remains there for ~ 150 000 
cycles. 




6 8 
MC cycles / 10 6 

FIG. 10. Qa during a canonical MC run at T = 0.19 efc" 1 . 
The label shows the short section of the simulation where the 
fee funnel is entered. 



The order parameter Q4 does not allow us to obtain the 
free energy barrier between the icosahedral and liquid- 
like states. However, it must be considerably smaller 
than the barrier for passage from the liquid-like state to 
the fee funnel, as dynamical coexistence of the two states 
is seen over a wide range of temperature. For LJ55 E c is 
able to detect a free energv-barrier between the Mackay 
icosahedron and the liquid,E3 but for LJ38 p(E c ) is uni- 
modal. However, p(Qe) has distinct maxima correspond- 
ing to icosahedral and liquid-like states in the region of 
the melting transition (Figure [ill) ; the separate max- 
ima disappear above T = 0.19efc _1 . At T = 0.18 efc" 1 
these maxima give rise to a free energy barrier of 0.93 fcT 
for passing from the liquid-like state into the icosahedral 
funnel. This compares to a value of 8.08 fcT for passing 
into the fee funnel. This difference is consistent with the 
picture of the energy landscape obtained from the dis- 
connectivity tree that showed the fee funnel to be much 
narrower (and therefore less accessible); it probably re- 
sults from the greater structural similarity between the 
icosahedral and liquid- like states. This result clearly ex- 
plains why the LJ38 cluster is much more likely to enter 
the icosahedral funnel on relaxation down the PES. 



IV. CONCLUSIONS 

The disconnectivity tree and free energy profiles that 
we have computed in this paper provide an integrated 
picture of the energy landscape of LJ 38 . The PES has two 
funnels. The fee funnel is deep and narrow and termi- 
nates at the truncated octahedral global minimum. The 
icosahedral funnel is much wider and has a flatter bot- 
tom. The large energy barrier between the two funnels 
gives rise to a large free energy barrier that causes the 



cluster to be trapped in one of the funnels when the tem- 
perature is below the melting point. Furthermore, the 
difference in the width of the funnels leads to a much 
lower free energy barrier for passage from the liquid- 
like state into the icosahedral funnel. This explains why 
the cluster preferentially enters the icosahedral funnel on 
cooling and why global optimization of this cluster is dif- 
ficult using those methods that do not transform the en- 
ergy landscape. 



liquid-like 




FIG. 11. Probability distribution of Qa during a MC run 
at T= 0.18 efc" 1 . 



LJ38 has a paradigmatic double-funnel energy land- 
scape and so our understanding of this cluster can help 
provide insight into other systems for which a number 
of funnels may be present on the PES. There is a clear 
relationship between multiple funnels and trapping, and 
we saw for LJ38 that the free energy barriers between 
low-energy structures in different funnels increase rela- 
tive to the thermal energy as the temperature decreases. 
In glass-forming systems we expect that the glass transi- 
tion is the result of similar increases in free energy bar- 
riers between different low-energy amorphous structures. 
It would therefore be very interesting to use disconnec- 
tivity trees to probe the multiple-funnel structure of the 
energy landscapes of model glasses. As water is a 'strong' 
liquid the disconnectivity tree that has been obtained far 
(H20)2o may provide a clue to what one might expecto 

The trapping that results from multiple funnels on the 
energy landscape also clearly shows why it.is important 
that proteins have a single accessible funnel.0 Trapping in 
a non-native funnel would prevent the protein from per- 
forming its biological function. Our results for LJ38 also 
suggest the possibility that the native state might not 
always correspond to the free energy global minimum. 
There could be some cases where the funnel leading down 
to the global minimum is so narrow that it is inaccessible 
on most time scales.Q The folding of the protein would 
then be associated with trapping in a more accessible fun- 
nel, and the functional lifetime of the protein would de- 
pend on the rate of escape from this funnel to the global 
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minimum. Indeed, there are some proteins from the ser- 
pin family of protease inhibitors for which this seems to 
occur E2 For example, the protein plasminogen activator 
inhibitor-1 first folds to the active state, but then on a 
time scale of hours .ran .spontaneously transform to an 
inactive latent form.E3H However, this scenario is likely 
to be the exception rather than the rule. 

We found that the ease of reaching the global mini- 
mum on a multiple- funnel energy landscape is related to 
the free energy barriers for entering the different funnels 
from the liquid, which in turn are related to the width 
of the funnels. Our results, therefore, add weight to the 
intuitive picture that glass-formers have a narrow funnel 
leading down to the crystal. These features of the energy 
landscape are intimately related to the structure. For ex- 
ample, in LJ38 the greater accessibility of the liquid-like 
state can also be explained in terms of the greater struc- 
tural similarity of the liquid to the icosahedral structures 
(both have some polytetrahedral character) than to the 
fee structures. Similarly, Straley found that crystalliza- 
tion to a bulk close-packed structure in flat space is much 
more difficult than to the completely polytetrahedral 
{3, 3, 5} polytope in a positively-curved space (the three- 
dimensional hypersurface of a four-dimensipnal sphere) 
where this structure is the global minimumEII 
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